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ABSTRACT 


This study proposes a new contribution to solve the problem of automatic 
landmarks detection in three-dimensional cephalometry. 3D images obtained 
from CBCT (cone beam computed tomography) equipment were used for 
automatic identification of twelve landmarks. The proposed method is based 
on a local geometry and intensity criteria of skull structures. After the step of 
preprocessing and binarization, the algorithm segments the skull into three 
structures using the geometry information of nasal cavity and intensity 


Keywords: information of the teeth. Each targeted landmark was detected using local 
3D cephalometry geometrical information of the volume of interest containing this landmark. 
The ICC and confidence interval (95% CI) for each direction were 0, 91 
Automatic algorithm (0.75 to 0.96) for x- direction; 0.92 (0.83 to 0.97) for y-direction; 0.92 (0.79 
CBCT image to 0.97) for z-direction. The mean error of detection was calculated using the 
Landmarks detection Euclidian distance between the 3D coordinates of manually and 
Orthodontics automatically detected landmarks. The overall mean error of the algorithm 
was 2.76 mm with a standard deviation of 1.43 mm. Our proposed approach 
for automatic landmark identification in 3D cephalometric was capable of 
detecting 12 landmarks on 3D CBCT images which can be facilitate the use 

of 3D cephalometry to orthodontists. 
This is an open access article under the CC BY-SA license. 
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1. INTRODUCTION 

Cephalometric analysis is the measurement of the relationship between reliable landmarks on the 
human skull. It is well-known for many applications, such as orthodontics and maxillofacial surgery. It was 
introduced for the first time in 1931 by taking measurements for different facial units using radiographic 
image [1, 2]. Cephalogram is widely used by orthodontics to perform manual cephalometric analyses based 
on manually identified landmarks [3-5]. The manual identification is tedious, time-consuming and needs high 
expertise for landmark localization. Thereafter several studies were proposed different methods for automatic 
landmark detection using radiographic images [6-9]. 

The main drawback of the said methods is that rendering a 2d representation of a 3d structure. 
CBCT (cone beam computed tomography) is the recent imaging technology with low radiation and cost 
compared to CT imaging. CBCT machines beginning to make a large attention in many areas of dental 
practice, such as implant placement and orthodontic diagnosis [10-12]. Gribel et al. [ 13] have demonstrated 
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that measurements in 3D CBCT are more accurate than lateral cephalogram. With CBCT we can get a precise 
representation of the anatomical structure; a study conducted by Swennen and Schutyserb [14] presented the 
advantages and disadvantages of CBCT and CT in 3D cephalometry. In 2006 Olszewski et al. [15] 
demonstrated that 3D analysis gave the same results and adequate diagnoses as 2D analysis using the same skull 
while Adam et al. [16] showed that using a 3D method was 4-5 times more accurate than the 2D approach. 

Nowadays automatic detection of hard landmarks using CBCT images presents the most challenge 
in 3D cephalometry. To solve this problem, recent studies have proposed various methods of locating 
landmarks automatically. In literature there are four main approaches adopted for automatic landmark 
localization: machine-learning [17], atlas-based image registration [18-20], model-based [21, 22] and 
knowledge-based [23, 24]. In 2014 Makram and Kamel [19] proposed a method for 3D landmark detection. 
Their method based on a reeb graph technique and mesh registration. In the same year, Shahidi et al. [20] 
proposed designed software for automatic landmark localization. They combined feature-based and intensity- 
based in image registration technique. Recently Codari et al. [18] proposed an approach for semi-automatic 
landmark detection in 3d cephalometry using automatic segmentation followed by an intensity-based 
registration technique. In 2015 an automatic algorithm using anatomical definition of landmarks was 
proposed by Gupta et al. [23]. 

However, this algorithm uses a static vectors for VOI extraction, which doesn’t take into account the 
anatomical variability and cases of deformed structures. A technical report published by Neelapu et al. [24] 
propose an automatic algorithm based on symmetrical features of skull to dynamically extract the region for 
each group of landmarks. In this work 11 landmarks are located in the midsagittal plan. Therefore, the 
detection of these landmarks using the proposed algorithm becomes difficult when the position of some 
landmarks do not aligned with the midsagittal plan. In study we aim to achieve some robustness against 
deviation of some landmark from the midsagittal plan, due to deformation and anatomical variability between 
structures. We propose an algorithm that automatically divides the skull into three parts: mandible, maxillae 
and the upper part containing orbit and nasal bone. The midsagittal plan of each structure was extracted 
followed by the landmark detection process of targeted landmarks the new method for landmark detection in 
3D cephalometry using geometric and intensity-based criteria of maxillo-facial structures and teeth. 


2. RESEARCH METHOD 
2.1. Dataset 

20 skulls of Moroccan people were used in this study. The scan of each patient was obtained by a 
CBCT scanner for other diagnostic reasons independent of this study. Each 3D dataset was converted into 
DICOM (digital imaging and communication in medicine) format with 609 2D slices (512*512 of spatial 
resolution) and voxel resolution of 0.3 mm. 


2.2. Gold standard generation 

Manual identification of each landmark was performed by three orthodontists. The identification 
was done using MPR (multi planar reconstruction) views (axial, sagittal and coronal plan). The mean 
coordinates for each landmark was used as gold standard to validate our algorithm. The location of each 
landmark was reported in the form of 3D Cartesian coordinates. 


2.3. Landmarks definitions 


Table 1 presents the definitions of 12 landmarks under consideration in this paper as defined by 
Swennen et al. [25]. 


Table 1. Anatomical landmarks definitions 


No Landmark Definition 

1 Menton: Men Menton is the most inferior midpoint of the chin on the outline of the mandibular symphysis. 

2 Gnathion: Gn The most anterior and inferior point on the contour of the mandibular symphysis. 

3 Pogonion: Pog Pogonion is the most anterior midpoint of the chin on the outline of the mandibular symphysis. 

4 B-Point: B B-Point is the point of maximum concavity in the mid- line of the alveolar process of the mandible. 

5 Anterior Nasal Anterior Nasal Spine is the most anterior midpoint of the anterior nasal spine of the maxilla. 
Spine: ANS 

6 A-Point: A A-Point is the point of maximum concavity in the mid- line of the alveolar process of the maxilla. 

ih Pronasal The most anterior point in nasal bone 

8 Nasion: N Nasion is the midpoint of the frontonasal suture. 

9-10 Orbitale: OrR-OrL Orbitale (Or) is the most inferior point of each infra-orbital rim. 

11-12  Frontozygomatic: Frontozygomatic (Fz) is the most medial and anterior point of each frontozygomatic suture at the 
FzR-Fz] level of the lateral orbital rim 


Bulletin of Electr Eng & Inf, Vol. 10, No. 2, April 2021: 707 -715 


Bulletin of Electr Eng & Inf ISSN: 2302-9285 o 709 


2.4. Overview of the automatic algorithm 

After the image preprocessing and segmentation step, another algorithm is used for the automatic 
separation of the skull into three main anatomical parts. Anatomical and intensity information from the nasal 
cavity and teeth is used to automate the entire process. The three parts are represented by the maxilla, the 
mandibular and the upper part of the skull. Figure 1 presents an overview of the necessary steps of the 
proposed algorithm. 
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Figure 1. Overview of the proposed algorithm 


2.5. Preprocessing and segmentation 

In this step the input volume was smoothed using a Gaussian filter. The size of this filter was set to 3 
in order to limit the blurring effect. Here is the Gaussian function for 2D image. After volume smoothing, a 
thresholding segmentation method was used to segment out the soft tissue in each volume. 


2 2 
X“ Fyr 


5 
26° 





G(x, y) = 





e 
2m0? 


2.6. Automatic partitioning of maxilla-facial structures 

In this step we aim at developing an automatic algorithm for partitioning the skull into three parts. 
Geometry and intensity are two main criteria used by the proposed algorithm. The anatomical and 
geometrical knowledge of nasal cavity were used to divide the skull into lower and upper part, while the 
intensity information of teeth was used to separate the lower part into maxillae and mandible partitions. 
Figure 2 shows an example of a maxillofacial structure automatically divided into three parts by the proposed 
algorithm. Figure 3 present an example of the three plans used to describe our proposed method for automatic 
landmark detection. Figure 2 shows a 3D surface rendering using the marching cube in a three dimensional 
coordinates in Cartesian system. In the following, we explain the three main steps of the proposed algorithm 
to separate the three structures. 





Figure 2. The estimated midsagittal plan with the three partitions of the skull 
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Figure 3. MPR (multi-planar reconstruction) 


2.6.1. Automatic detection of anterior profile 


The steps used can be presented as: 
The algorithm starts by a global binarization using a single threshold equal to 700 grey levels for the 
whole skull. The 700 value was chosen based on experimental test using the 20 subject under study. 
Detection all voxels with minimum x-coordinate in each axial slice. 
The sagittal plan represented by mean of y-coordinate values of all detected voxels, was estimated as 
midsagittal plan. 
Extraction of anterior profile of the estimated midsagittal plan by selecting the voxels with minimum x- 
coordinate in each axial slice. The Figure 4 represents the estimated anterior profile. 
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Figure 4. Anterior profile of estimated midsagittal plan figure 


2.6.2. Nasal cavity detection algorithm 


The localization of center of nasal cavity is done by analyzing the anterior profile of midsagittal plan 


based on anatomical and geometrical information of anterior profile skull. 


The analysis starts by determining the location of maximum rate change. The first derivative of the 
function represented the anterior profile allows finding these rates , see Figure 5, which are characterized 
by a value greater than 100 or less than -100 (empirical values). To get better results of differentiation all 
images have been pre-processed using smoothing techniques. For each point with maximum rate change, 
find the point with x-coordinate value less than the x-coordinate value of this point, then calculate the size 
of the segment between these two points. 

The segment with maximum size tends to represent the nasal cavity. Figure 5 shows the result segment in 
green that represent the location of nasal cavity. 

Find the point in the middle of the reported segment as shown in Figure 5 that tends to estimate the center 
of nasal cavity. 
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Figure 5. Gradient of anterior profile 


2.6.3. Maxillary and mandibular alveolar estimation algorithm 

— After Anterior profile detection, find the point with maximum intensity in the anterior profile which 
represents the seed of teeth region. 

— Select all voxel in lower and upper side that represent a difference of 100 gray levels with seed voxel. 

— Select the median voxel in the detected region. 

— Median voxel plus 40 voxel (empirical value) estimate maxillary alveolar process while Median minus 40 
voxel estimates mandibular alveolar process. 


2.7. Contour extraction 
The contour in each partition is computed by finding the minimum x-coordinate values in each axial 
slice. The technique of contour extraction has been described in previous work by Gupta et al. [23]. 


2.8. Landmark detection 
After contour extraction, landmark detection is based on the definition given in Table 1. 


2.8.1. Landmark Menton (Me), Gnathion (Gn), Pogonion (Pog) and B-point (B). 

These landmarks are located in the midline of the symphysis menti. The algorithm uses the 
symphysis menti contour to locate these landmarks. First a line segment was formed using first and last point 
on of extracted contour, then the distances between the line segment and each point in the contour was 
calculated. Based on these distances and the anatomical definition of the landmarks in Table 1, the four 
landmarks were detected, as shown in Figure 6. The maximum distance between the segment line and the 
contour under line is considered as the Pogonion landmark and the maximum distance from contour above 
line segment is detected as B-point. The first point on the contour is the Menton landmark and the point in 
the middle of Pogonion and Menton is considered as Gnathion landmark. 
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Figure 6. Extracted contour in black, line segment in green and detected landmarks in red based on distances 
from line segment 
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2.8.2. Anterior nasal spine: ANS and A-point: A 

After the contour detection step, a fully automatic algorithm is designed based on local anatomical 
property of these landmarks. The algorithm start by finding the horizontal line segment in the contour and 
bounded by two distinct end point and characterized by: have a size greater than 10, and x- coordinate of 
second point less than x-coordinate of first point. Figure 7 shows contour, line segment and detected 
landmarks. Ans is detected as the second end point and A-point considered as the point with largest 
perpendicular distances between line segment and contour. 
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Figure 7. Contour and line segment in green for Ans and A-point (red point) landmarks detection 


2.8.3 Pronasal landmark 

This point is the most anterior point in the nasal bone. By crossing the upper partition of 3D volume 
from anterior to posterior direction, nasal bone was the first structure that had founded. The first point of 
intersection between nasal bone and the coronal plane is considered as the Pronasal landmark as Figure 8(a) 
shows. 


2.8.4. Nasion: N 

The same technique used for detection Ans and A-point is used for Nasion landmark. Horizontal line 
segments that have size of 10 and x-coordinate of second point less or equal to x- coordinate of first point. 
The middle point of line segment was detected as Nasion. Figure 8(b) shows contour, line segment and 
Nasion landmark. 
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(a) (a) 
Figure 8. (a) Poronasal showed in sagittal plan with and without soft tissue, (b) contour and line segment for 


Nasion (red point) detection 


2.8.5. Orbitale: OrR-OrL and Frontozygomatique: FzR-Fzl 
The algorithm proposed for Orbital landmark identification take Pronasal landmark as input to start 
the detection process. The steps of this algorithm are presented as follows: 
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— Step 1: Using the axial slice of Pronasal landmark, the algorithm finds the y-coordinate of orbit center 
using the y-coordinate of most lateral voxel and y-coordinate of Pronasal landmark. 

— Step 2: using the orbit center, a region of interest automatically selected by the proposed algorithm. 

— Step 3: determination of Frontozygomatic landmark starting from lateral orbital rim where the lateral 
orbital rim is the most anterior point as shown in Figure 9. 

— Step 4: Starting from the orbit center, the algorithm traverses the 3D volume of skull in the lower 
direction until it found a closed contour that represent the orbital landmark. Figure 10 an example of axial 
slices used for the detection process. 





Figure 9. Step 3-estimated Figure 10. Step 4 axial slices used to detect orbital landmark 
Frontozygomatic landmark 
in coronal slice 


2.9. Data analysis 

In this study, the proposed method was testes on 20 CBCT images. ICC and their 95% confident 
intervals were used to assess interobserver reliability for each landmark in the x-, y- and z-direction. 
Thestatistical test done by IBM SPSS statistical package version 23 (SPSS Inc., Chicago, IL) based on a 
mean-rating (k=3), absolute-agreement, 2-way random-effects model. 


3. RESULTS 

The ICC and confidence interval (95% CI) for each direction were 0.91 (0.75 to 0.96) for x- 
direction; 0.92 (0.83 to 0.97) for y-direction; 0.92 (0.79 to 0.97) for z-direction. 3D Euclidian distance 
formula were used to calculate the detection error between manual and automatic landmark identification. 
The overall mean error of landmark detection was 2.76 mm. 





Distance= f (x1 — x2)}* + (yl — y2} + (21 — 22)" 


Table 2 shows the mean error between automatic and manual landmark detection. The overall mean error 
was equal to 2.76(+1.40) mm, with 1.34 (40.77) mm as minimum error and 4.59 (+2.35) mm as maximum 
error. 


Table 2. Mean error and standard deviation for each landmark 


N Mean error (mm) STD 
1 1.97 1.30 
2 2.76 1.75 
3 2.93 1.45 
4 4.59 2.35 
5 4.01 1.52 
6 1.34 0.77 
7 1.80 0.99 
8 3.86 2.28 
9 3.03 0.99 
10 3.30 1.53 
11 1.57 1.08 
12 1.98 0.77 
Mean 2.76 1.40 
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4. DISCUSSION 

This paper has proposed a new algorithm for automatic landmark detection on 3D maxillo-facial 
structures. The authors have focused on twelve landmarks: Menton, Gnathion, Pogonion, B-point, Ans, A- 
point, Orbitale left and right, Frontozygomatique left and right, Pronasal and Nasion. Twenty CBCT images 
were used to validate our algorithm. The proposed method was capable of locating 12 landmarks in all 
subjects under study. Euclidean distance was used to calculate the mean error between manual and automatic 
landmark detection. Analysis of inter-observer reliability was done using ICC index reliability following the 
guideline proposed by T. Koo and M. Li [26]. The ICC and 95% confidence intervals (95% CI) for inter- 
observer reliability was equal to 0.92 (0.79 to 0.96). It is interesting to note that the error for bilateral 
landmarks is significant compared to the error for midsagittal landmarks. This difference may explain by the 
ambiguity in location of these landmarks by an operator using multi-planar reconstruction view only. B. 
Hassan et al. [27] find that performing manual landmark detection using MPR (multiplanar reconstruction) 
view with combination of 3D surface models improve precision but increase time of detection. 

The errors reported by our method can be compared with the mean errors of state of the art [17-20, 
22, 23, 28] which ranging from 1.2 mm to 3.64 mm with a maximum error of 9.7 (+3.8) mm. Table 2 shows 
that 58.3 % of landmarks were detected with a mean error less than 3 mm, 83.3% of landmarks have a mean 
error less than 4 mm and 100% of landmarks identified with mean error less than 5 mm. Image denoising is 
very important setp to improve the quality of image [29], many studies doesn’t explain the technique used for 
noise reduction. The most of existing works based on model , training or atlas approach to develop automatic 
algorithm. These technique need more resources and data which limits the algorithms toward use in clinical 
practice. Neelapu et al. [24] state that atlas-based technique is prone to errors due to skeletal deformiries. In 
this study the proposed algorithm use only landmarks definition without any training step or existing model. 
In the best of our knowledge, there is one study conducted by Gupta et al. [23] that addressed the approach of 
knowledge-based for automatic landmark detection in 3D cephalometry. In the previous work the algorithm 
proposed is based on static predefined vector to define volume of interest while in our study the volume of 
interest was defined dynamically for each image. All the twenty images were successfully separated into 
three main structures by the algorithm using a combination of geometry and intensity criteria. Automatic 
partitioning of the skull lead to solve tow problems that exist in the previous works. The first problem is 
landmarks detection of midsagittal structures in case of skeletal deformity and the second problem is human 
inter-individual anatomical variability. In this study a midsagittal plan was locally extracted for each 
structure and the identification of the volume of interest containing the target landmark was extracted 
dynamically by incorporating geometrical local information without any preset static information. We 
believe that using only anatomical geometry information of maxillo-facial structures can reduce the 
computational complexity of the algorithm compared to other approaches such as model-based and atlas- 
based. 


5. CONCLUSION 

In this study we have proposed an automated method for automatic landmark detection in 3D 
cephalometry. We hope that our study will valuable for solving the difficulty of 3D cephalometric analysis 
using CBCT images. One limitation of CBCT image in cephalometry is still the large field of view, which 
lead to higher radiation dose for the patient. Despite this, we believe that our method could help researchers 
aiming at developing automatic system for 3D cephalometry using CBCT modality. 
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